

/* Procedure that gives back the fractional index for a new value
   in a table. This is a sort of linear interpolation */

using namespace std;
extern vector<double> tab_mass_axis;
//double get_star_index(vector<double> table_input, double value_evaluate, 
double get_star_index(double value_evaluate, 
	long ntime_axis, double tinx_interp, double last){
 
//  double dum_mass_axis[nmass_axis]; 
//  for(long tt=0;tt<(long)nmass_axis;tt++){
 //     dum_mass_axis[tt]=table_input[tt*ntime_axis+(long)tinx_interp];
  //}


  //Input array is monot decreasing
  //Start from the tail, most likely to be there 
  long index, val3, val5;
  val5=(long) ceil(last);
  //find your location
  val3=(val5+1)*ntime_axis+(long)tinx_interp;
//cout<<"calc ind\t"<<last<<"\t"<<ntime_axis<<"\t"<<tinx_interp<<endl;
//  if((long)val3 > (long)tab_mass_axis.size()){
//	cout<<val3<<"\t"<<tab_mass_axis.size()<<endl;
//	 exit(0);
 // }
  for(index=val5;index>=0;index--){
    val3-=ntime_axis;
//cout<<val3<<"\t"<<tab_mass_axis.size()<<endl;
//    if(value_evaluate > 119) cout<<tab_mass_axis[val3]<<"\t"<<val3<<"\t"<<index<<endl;
    if(tab_mass_axis[val3] > value_evaluate) break;
  }
 
  
  //get the fract index
  double fract_index;
if (index > 0){
  double val1, val2;
  val3=(long) index*ntime_axis+(long) tinx_interp;
  val1=tab_mass_axis[val3];
  val2=tab_mass_axis[val3+ntime_axis];
  fract_index=(val1-value_evaluate)/
  	(val1-val2)+index;
  }else(fract_index=0.0);
// index=0;  //fix this kludge!!!! -RDS
//if (value_evaluate > 119.){
//  cout<<fract_index<<"\t"<<value_evaluate<<"\t"<<index<<"\t"<<tab_mass_axis[index]<<"\t"<<tab_mass_axis[index+1]<<endl;
//  }
//if(value_evaluate > 119) cout<<"exit gsi ok"<<endl;
 
  return fract_index;
  
}
